Report

library(plotly)
library(Hmisc)
library(WVPlots)
library(plyr)
source("tools/DataSplitTools.R")
source("tools/GeneralizedNadarayaWatson.R")
source("tools/CommonTools.R")
source("tools/CrossValTools.R")

Аналіз даних

Зчитаємо дані:

df <- read.csv2(file = "data/ZNOandVoating/input.csv", 
                header = FALSE, sep = ";", dec = ",")
names(df) <- c("ukr","math", "pro_ukr", "radical", "opposition", "small", "not_voted")

head(df, 3)

##   ukr math  pro_ukr  radical opposition     small not_voted
## 1 112    0 19.37149 3.831828  10.466742 13.519935     52.81
## 2 133    0 41.61113 8.886996   1.342683  8.369190     39.79
## 3 149  125 21.52264 4.991798  11.615622  9.729938     52.14

ЗНО з математики та української мови

Виведемо загальні характеристики даних ЗНО з математики та української мови:

describe(df[, c("math", "ukr")])

## df[, c("math", "ukr")] 
## 
##  2  Variables      268003  Observations
## ---------------------------------------------------------------------------
## math 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   268003        0       55    0.774    53.95    70.16        0        0 
##      .25      .50      .75      .90      .95 
##        0        0      125      157      172 
## 
## lowest :   0 100 104 107 111, highest: 196 197 198 199 200
## ---------------------------------------------------------------------------
## ukr 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   268003        0       83    0.997    125.3    58.57        0        0 
##      .25      .50      .75      .90      .95 
##      111      137      165      184      189 
## 
## lowest :   0.0 100.0 101.0 103.0 104.0, highest: 197.0 198.0 199.0 199.5 200.0
## ---------------------------------------------------------------------------

Із таблиці вище видно, що в записах нема пропущених значень. Судячи із значень квантилів, половина значень ЗНО з математики є нулями. Варто відмітити, що значно менша частка значень 0 у ЗНО з української мови.

Подивимося на розподіли:

p1 <- plot_ly(x = df[, "math"], type = "histogram") %>% layout(title = "EIT: Mathematics")
p2 <- plot_ly(x = df[, "ukr"], type = "histogram") %>% layout(title = "EIT: Ukranian language")

plotly::subplot(p1, p2)

remove(p1); remove(p2)

Подивимося на співвідношення людей, котрі отримали 0 по математиці та українській одночасно

print("Не подолали поріг ЗНО відносно всіх абітурієнтів: " %&% as.character(round(nrow(df[(df$math == 0 & df$ukr == 0), ]) * 100 / nrow(df), 2)) %&% "%")

## [1] "Не подолали поріг ЗНО відносно всіх абітурієнтів: 12.66%"

print("Не подолали поріг ЗНО відносно тих, хто не склав математику: " %&% as.character(round(nrow(df[(df$math == 0 & df$ukr == 0), ]) * 100 /
                                                                                     nrow(df[(df$math == 0), ]), 2)) %&% "%")

## [1] "Не подолали поріг ЗНО відносно тих, хто не склав математику: 20.8%"

print("Не подолали поріг відносно тих, хто не склав українську: " %&% as.character(round(nrow(df[(df$math == 0 & df$ukr == 0), ]) * 100 / 
                                                                                       nrow(df[(df$ukr == 0), ]), 2)) %&% "%")

## [1] "Не подолали поріг відносно тих, хто не склав українську: 91.04%"

Із останнього напрошуєтся висновок, що ті, хто не здали українську, скоріш за все не здали й математику.

Поглянемо на кореляції між результатами ЗНО, політичними силами:

c <- cor(df, method = "pearson")
plot_ly(x = rownames(c), y = colnames(c), z = c, colors = "Greys", type = "heatmap") %>% layout(title = "Pearson Correlation")

remove(c)

Відмітимо майже 0 кореляцію між результатами ЗНО з української мови та про-український політичний напрямок та слабку кореляцію між результатми математики та української.

Видалимо стрічки, що відповідають спостереженням із 0 хочаб в одному ЗНО та подивимося на розподіл.

df <- df[(df$math != 0) & (df$ukr != 0) ,]
p1 <- plot_ly(x = df[, "math"], type = "histogram", name = "math")
p2 <- plot_ly(x = df[, "ukr"], type = "histogram", name = "ukr")

plotly::subplot(p1, p2) %>% layout(title = "Distribution of the passed EIT",
  annotations = list(list(x = 0.2, y = 1.05, text = "math", showarrow = F, xref = 'paper', yref = 'paper'),
                     list(x = 0.8, y = 1.05, text = "ukr", showarrow = F, xref = 'paper', yref = 'paper')))

remove(p1); remove(p2);

Відмітимо, що наразі дуже відсіялася частка людей, котрі погано здали українську мову.

Гістограма розсіювання, в свою чергу, виглядає наступним чином:

plot_ly(data = df, x = ~math, y = ~ukr, mode = "markers", type = "scatter") %>% layout(title = "Passed EIT results")

Варто відмітити, що для абітурієнтів, що склали математику краще ніж 160 балів, помітна тенденція, як правило, кращого складання ЗНО з математики.

Подивимося на розподіли по кожній політичній силі окремо (якщо конкретно обрана політична сила має для даного спостереження найбільший вплив, то дане спостереження відносимо до ціє політичної сили)

• математика

W <- as.data.frame(df[, -c(1, 2)])
W[, "max_value"] <- apply(W, 1, max)

all_plots <- list()
annotations_list <- list()
names_ <- names(W[, -ncol(W)])

for (i in seq_along(names_)) {
  name <- names_[i]
  d <- df[W[, name] == W["max_value"], ]
  if (nrow(d) != 0) {
    all_plots[[name]] <- plot_ly(x = d[, "math"], type = "histogram", name = name)
    annotations_list[[i]] <- list(x = 0.15*i, y = 1.005, text = name, showarrow = F, xref = 'paper', yref = 'paper')
  }
}

plotly::subplot(all_plots) %>% 
  layout(title = "Mathematics EIT according to political strength", annotations = annotations_list)

remove(d); remove(all_plots); remove(annotations_list)

• українська мова

all_plots <- list()
annotations_list <- list()

for (i in seq_along(names_)) {
  name <- names_[i]
  d <- df[W[, name] == W["max_value"], ]
  if (nrow(d) != 0) {
    all_plots[[name]] <- plot_ly(x = d[, "ukr"], type = "histogram", name = name) 
    annotations_list[[i]] <- list(x = 0.15*i, y = 1.005, text = name, showarrow = F, xref = 'paper', yref = 'paper')
  }
  
}

plotly::subplot(all_plots[["pro_ukr"]], all_plots[["not_voted"]]) %>% 
  layout(title = "Ukr EIT according to political strength", annotations = annotations_list)

remove(d); remove(all_plots); remove(annotations_list)

Відмітимо, що розподіл результатів ЗНО з математики не на стільки різниться від розподілу результатів з української. Також варто відмітити, значно меншу кількість абітурієнтів, котрі склали ЗНО з укр. мови в межах 100 з областей, що імовірніше підтримали про-українські позиції.

all_plots <- list()
for (name in names(W[, -ncol(W)])) {
  d <- df[W[, name] == W["max_value"], ]
  if (nrow(d) != 0) {WVPlots::ScatterHist(d, "math", "ukr", name, smoothmethod = "none", estimate_sig = FALSE,  minimal_labels = TRUE)}
}

Із останньої пари зображень здається, що нижній поріг результатів української (для абітурієнтів із математикою > 160 балів) вищий у випадку про-українських регіонів.

Результати голосування

Виведемо загальні характеристики кожної політичної сили:

df <- read.csv2(file = "data/ZNOandVoating/input.csv", 
                header = FALSE, sep = ";", dec = ",")
names(df) <- c("ukr","math", "pro_ukr", "radical", "opposition", "small", "not_voted")
describe(df[, -c(1, 2)])

## df[, -c(1, 2)] 
## 
##  5  Variables      268003  Observations
## ---------------------------------------------------------------------------
## pro_ukr 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   268003        0       22    0.994    34.13    14.61    15.71    16.23 
##      .25      .50      .75      .90      .95 
##    21.52    33.81    43.01    53.87    53.87 
## 
## lowest : 15.71244 16.22691 19.37149 20.13312 21.00261
## highest: 43.00824 44.77893 48.82355 50.41795 53.87200
## ---------------------------------------------------------------------------
## radical 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   268003        0       22    0.994    6.819    2.525    3.047    3.832 
##      .25      .50      .75      .90      .95 
##    4.992    8.078    8.477    8.912   10.428 
## 
## lowest :  3.046992  3.831828  3.834072  4.523960  4.619912
## highest:  8.886996  8.911950 10.084956 10.427880 11.555155
## ---------------------------------------------------------------------------
## opposition 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   268003        0       22    0.994    4.482    4.851   0.4165   0.4970 
##      .25      .50      .75      .90      .95 
##   1.0428   1.8701   7.1334  11.6156  14.5749 
## 
## lowest :  0.344142  0.416508  0.497000  0.953295  1.042825
## highest:  6.796640  7.133360 10.466742 11.615622 14.574912
## ---------------------------------------------------------------------------
## small 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   268003        0       22    0.994    9.201    2.142    6.779    7.189 
##      .25      .50      .75      .90      .95 
##    7.553    8.845    9.966   13.113   13.113 
## 
## lowest :  5.856787  6.778902  7.189344  7.360584  7.496422
## highest: 10.890088 11.198572 11.346280 13.112736 13.519935
## ---------------------------------------------------------------------------
## not_voted 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##   268003        0       22    0.994    45.37    10.68    30.00    30.00 
##      .25      .50      .75      .90      .95 
##    39.79    45.36    52.81    57.20    60.48 
## 
## lowest : 30.00 31.72 35.15 36.27 39.79, highest: 54.68 55.32 57.20 58.64 60.48
## ---------------------------------------------------------------------------

Одразу видно, що 22 унікальні значення в кожному стовпчику відповідають 22 областям України. Варто відмітити, що відсутні області, із переважними радикальними або опозиційними силами. Також відмітимо що немає регіонів, у яких голоса віддані виключно за малі політичні партії.

Розподіл абітурієнтів по областях:

ab <- plyr::ddply(.data = df, .variables = .(pro_ukr, radical, opposition, small, not_voted), .fun = nrow)
plot_ly(ab, y = ~V1, type = "bar") %>% 
  layout(title = "Abiturients distribution according regions", 
         xaxis = list(title = "Region ID"), yaxis = list(title = ""))

Варто відмітити, що два регіони надто відрізняються чисельністю. Через таку зміщенність даних можуть бути зміщені результати.

Із розподілу відсотків підтримки кожної політичної сили відповідно до регіону видно, що більшість регіонів є політично неактивними. Також варто відмітити, що розподіли голосів малих партії, у більшості випадків, збігається з розподілом голосів радикальної партії.

plot_ly(y = ab[, "pro_ukr"], type = "bar", name = "pro_ukr") %>%
  add_trace(y = ab[, "not_voted"], name = "not_voted") %>%
  add_trace(y = ab[, "radical"], name = "radical") %>%
  add_trace(y = ab[, "opposition"], name = "opposition") %>%
  add_trace(y = ab[, "small"], name = "small")

remove(ab)

Висновок аналізу даних

На основі порівняння розподілів можна зробити висновок про те, що дані можуть бути описані за допомогою моделі суміші зі змінними концентраціями. Проте, зважаючи на кореляції та відсоткові співвідношення, можливо, модель суміші зі змінними концентраціями, буде не надто еффективною.

Розбиття

Розіб’ємо вибірку на тренувальну та тестову частини у відсотковому співвідношенні 80/20%. Далі, тренувальну частину розіб’ємо на 5-частин для проведення 5-ти фолдової кроссвалідації з метою визначення оптимального для кожної компоненти параметра згладжування h. Узагалі кажучи, можливо, варто було б врахувати незбалансованість даних під час розбиття, проте ми цього не робили.

df <- read.csv2(file = "data/ZNOandVoating/input.csv", 
                header = FALSE, sep = ";", dec = ",")
names(df) <- c("ukr","math", "pro_ukr", "radical", "opposition", "small", "not_voted")

df[, -c(1, 2)] <- df[, -c(1, 2)]/100
df <- df[(df$math != 0) & (df$ukr != 0) ,]

train_test_list <- train_test_split(df, ratio = 0.80)

train <- train_test_list[["train"]]
test <- train_test_list[["test"]]

cv_split <- cross_validation_split(train)
remove(df); remove(train_test_list)

Узагальнені оцінки Надарая-Ватсона

Визначення параметрів згладжування

Скористаємося 5-ти фолдовою кросс-валідацією для знаходження оптимальних параметрів згладжування для кожної моделі. Для визначення оптимального параметра використаємо МНК зважений модулем вагових коефіцієнтів:

\[MSE^{(m)} = \frac{1}{N} | \sum_{i=1}^{N}{a_{i:N}^{(m)}(Y_i - \hat{g}^{(m)}(X_i))^2} |\] Зауважимо, що модуль використовується для того щоб підчас підрахунків середнього зваженого MSE по всіх фолдах фіксованої компоненти, фіксованого h, не отримати 0 (який може бути хибно протрактований).

Отже, визначення оптимальних параметрів h відбувається за схемою описаною нижче.

• Зафіксуємо h, порахуємо зважений МНК для кожного фолда та знайдемо середнє та дисперсію по 5-ти фолдам для кожної компоненти. (Так прорахуємо для кожного фіксованого h.)

• Відкинемо 25% значень середніх, які відповідають найбльшим значенням дисперсії. Таким чином залишаємо лише ті результати, на яких моделі показали більш стабільну поведінку.

• Серед значень середніх результатів моделей визначаємо h, які відповідають найменшим середнім зваженим МНК (окремо для кожної компоненти).

Результат середніх та диспресії збережемно в змінній \(GNW_cv_results\), а оптимальні \(h\), серед запропонованих для перебору, у змінну \(h_opt\).

h_range <- c(seq(0.1, 1, 0.1), seq(1.5, 5, 0.5))

GNW_cv_results <- GNWcv_across_h(h_range = h_range, 
                                 cv_df_split = cv_split, 
                                 X_colname = "math", Y_colname = "ukr", 
                                 W_colname = c("pro_ukr", "radical", "opposition", "small", "not_voted"), 
                                 use_parallel = TRUE)

## [1] "h = 0.1"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] "h = 0.2"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] "h = 0.3"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] "h = 0.4"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] "h = 0.5"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] "h = 0.6"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] "h = 0.7"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] "h = 0.8"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] "h = 0.9"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] "h = 1"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] "h = 1.5"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] "h = 2"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] "h = 2.5"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] "h = 3"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] "h = 3.5"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] "h = 4"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] "h = 4.5"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] "h = 5"
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5

h_opt <- optimal_h(GNW_cv_results)

remove(cv_split); remove(h_range)

Збережемо про всяк випадок результати середніх та дисперсій кросс-валідації.

write.csv(GNW_cv_results, file = "data/computation_results/5_cv_mean_std.csv", row.names = TRUE)

Порівняння результатів кросс-валідації та прогнозування на відкладеній вибірці

# GNW_cv_results <- read.csv("data/computation_results/5_cv_mean_std.csv",  row.names = NULL, header= TRUE)
# names_ <- GNW_cv_results[, "X"]
# GNW_cv_results <- as.matrix(GNW_cv_results[, -1], mode="numeric")
# rownames(GNW_cv_results) <- as.vector(names_, mode="character")
# GNW_cv_results
# 
# h_opt <- optimal_h(GNW_cv_results)
# h_opt

Порівняємо тепер результати кросс-валідації з результатми прогнозування на тестовій частині даних.

GNW <- GeneralisedNadarayaWatson$new()
GNW$train(X_train = as.numeric(train[, "math"]), 
          Y_train = as.numeric(train[, "ukr"]), 
          W_train = as.matrix(train[, -c(1, 2)]))
GNW$run_cluster()
prediction_with_acoeff <- GNW$predict_in_parallel(X_test = as.numeric(test[, "math"]), 
                                                  W_test = as.matrix(test[, -c(1, 2)]), h = h_opt)
GNW$stop_cluster()

prediction <- prediction_with_acoeff[["prediction"]]
A_coeff <- prediction_with_acoeff[["A_test"]] 

remove(GNW); remove(train)

Порахуємо зважену суму МНК для прогнозу зробленого на тестових даних.

WMSE <- weighted_MSE(Y_true = as.numeric(test[, "ukr"]), Y_predicted = prediction, A_coeff = A_coeff)

Відмітимо, що найменше значення зваженого МНК для моделей, що відповідають проукраїнській силі та тим, хто не прогоолосував. Нагадаємо, найбільша по модулю кореляція спостерігалася між результатами ЗНО з української мовои та проукр. силами й тими, хто не проголосував. Також відмітимо, що раніше ми змогли виділити лише регіони, де переважає проукр напрямок та ті, в яких переважає не голосування. Тобто не виявилось регіонів, у яких переважають голоси за радикалів, за опозицію або за малі партії.

Порівняємо з результатами крос-валідації з відповідними параметрами згладжування.

cv_mean_results <- sapply(1:length(h_opt), 
                          function(i){
                            GNW_cv_results[which("mean_" %&% as.character(h_opt[i]) == rownames(GNW_cv_results)), i]
                            })
cv_test_comparing <- rbind(cv_mean_results, WMSE)
rownames(cv_test_comparing) <- c("cross-validation", "test")
colnames(cv_test_comparing) <- 1:5
print(cv_test_comparing)

##                         1        2        3         4          5
## cross-validation 507.7665 503.2183 353.0815  1229.091   295.0415
## test             454.6054 401.0795 299.7595 39133.250 -2553.6710

Із останньої таблиці видно, що квадрати залишків першої моделі та останньої, порахованих для кросс-валідації та на відкладеній тестовій вибірці не надто різняця у порівнянні з іншими компонентами.

Це свідчить про більш-менш вдалий підбір параметрів згладжування для моделей, що відповідають проукраїнському настрою та байдужому.

Для 3 інших моделей результати невтішні. Останнє може бути зумовлене слабкою представленістю відповідних класів у вибірці. Можливо є сенс розглянути 3-х компонентну суміш: проукраїнську направленість, “не проукраїнську” (об’єднати результати голосуваня за малі партії, опозицію, радикалів) та “байдужих”.

Візуалізація результатів

Зобразмо спочатку загальну картину: гістограму розсіювання з усіма прознозами для всіх моделей:

plot_ly(data = test, x = ~math, y = ~ukr, mode = "markers", type = "scatter", name="EIT results") %>% 
  add_trace(x=sort(test[, "math"]), y = prediction[sort(test[, "math"], index.return = TRUE)$ix, 1], name = "pro_ukr",  mode = 'lines+markers') %>%
  add_trace(x=sort(test[, "math"]), y = prediction[sort(test[, "math"], index.return = TRUE)$ix, 2], name = "radical",  mode = 'lines+markers') %>%
  add_trace(x=sort(test[, "math"]), y = prediction[sort(test[, "math"], index.return = TRUE)$ix, 3], name = "opposition",  mode = 'lines+markers') %>%
  add_trace(x=sort(test[, "math"]), y = prediction[sort(test[, "math"], index.return = TRUE)$ix, 4], name = "small",  line=list(color="yellow"), marker=list(color="yellow", line=list(color="yellow")), mode = 'lines+markers') %>%
  add_trace(x=sort(test[, "math"]), y = prediction[sort(test[, "math"], index.return = TRUE)$ix, 5], name = "not_voted",  mode = 'lines+markers') %>%
  layout(title = "Passed EIT results (test dataset)")

Окремо зобразимо для абітурієнтів, які проживають в областях, що ймовірніше наслідують проукраїнську групу або не голосують:

Wtest <- as.data.frame(cbind(prediction, test[, -c(1, 2)]))
Wtest[, "max_value"] <- apply(Wtest[, -(1:5)], 1, max)

pro_ukr <- test[Wtest[, "pro_ukr"] == Wtest["max_value"], ]
p1 <- plot_ly(data = pro_ukr, x = ~math, y = ~ukr, mode = "markers", type = "scatter", name = "EIT results in pro_ukr regions") %>% 
  add_trace(x = test[rownames(pro_ukr), "math"], 
            y = Wtest[rownames(pro_ukr), 1], name = "pro_ukr model",  mode = 'markers')


not_voted <- test[Wtest[, "not_voted"] == Wtest["max_value"], ]

p2 <- plot_ly(data = not_voted, x = ~math, y = ~ukr, mode = "markers", type = "scatter", name = "EIT results in not voted regions") %>%
  add_trace(x = test[rownames(not_voted), "math"],
            y = Wtest[rownames(not_voted), 5], name = "not_voted model",  mode = 'markers')

plotly::subplot(p1, p2)

rm(list = setdiff(ls(), lsf.str()))

Відмітимо, що основні тренди моделі для проукраїнських регіонів збереглися і лише кілька прогнозів вийшли надто далеко за межі допустимих оцінок (0 - 200).